Functional Connectivity Fingerprinting

2020-06-25 23:55:22

Author: Cedric Huchuan Xia (email, github)

Affiliation: Penn Lifespan Informatics and Neuroimaging Center (PennLINC)

---

1. Setup Environment

require(reshape2)
require(sna)
require(tidyr)

project_path = "~/Documents/xia_gps/"
data_path = file.path(project_path,"data/flywheel_data/network_txt")

2. Setup Target and Database FC

atlases = c("HarvardOxford","desikanKilliany","aal116","schaefer200x7","power264","gordon333","glasser360","schaefer400x7")
subj_net_files_1 = list()
subj_net_files_2 = list()

for (atlas in atlases) {
  subjects = list.files(file.path(data_path))
  for (subj in subjects) {
    subj_files = list.files(file.path(data_path,subj))
    net_pattern_1 = paste0("*task-rest*single*",atlas,"*")
    net_pattern_2 = paste0("*task-rest*multi*",atlas,"*")
    
    file_path_1 = file.path(data_path, 
                            subj,subj_files[grep(glob2rx(net_pattern_1),subj_files)])
    file_path_2 = file.path(data_path, 
                            subj,subj_files[grep(glob2rx(net_pattern_2),subj_files)])
    
    if (length(file_path_1)>0) {subj_net_files_1[[atlas]][[subj]] = read.table(file_path_1) }
    if (length(file_path_2)>0) {subj_net_files_2[[atlas]][[subj]] = read.table(file_path_2) }
  }
}

Visutalize FC matrix

matrix_fc = matrix(NA, 264,264)
matrix_fc[lower.tri(matrix_fc, diag = F)] = subj_net_files_1$power264[[1]][,1]
matrix_fc = symmetrize(matrix_fc, rule = "lower")
levelplot(matrix_fc,scales=list(draw=FALSE),col.regions = rev(rainbow(80))[-c(1:5)], region =T, ylab.right = "Pearson correlation", main=list(label='Singleband FC'),xlab="",ylab="")

4. Create Subj Similarity Matrix

gps_wide_matrix = data.frame()

for (subj in arrange(subj_days,`Days Collected`)$IID){
  #if ((subj %in% subj_days$IID[which(subj_days$`Days Collected`<=10)]) == T) {
    for (part in 1:1){
      half_1 = gps_clean2_feature$subj_mat_1[[subj]][[part]]$cor
      half_2 = gps_clean2_feature$subj_mat_2[[subj]][[part]]$cor
      half_1_half_2 = c(half_1,half_2)
      gps_wide_matrix = rbind(gps_wide_matrix,half_1)
      gps_wide_matrix = rbind(gps_wide_matrix,half_2)
      }
    #}
}

gps_wide_matrix = t(gps_wide_matrix)
gps_corplot = rquery.cormat(gps_wide_matrix, type = "full",graph=FALSE)
gps_corplot$subj = arrange(subj_days,`Days Collected`)$IID

levelplot(gps_corplot$r,scales=list(draw=FALSE),col.regions = rev(rainbow(1000))[-c(1:20)], region =T, ylab.right = "Pearson correlation", main=list(label='GPS Feature Similarity'),xlab="",ylab="")

fc_corplot$schaefer400x7$plot

fc_corplot$schaefer200x7$plot

fc_corplot$power264$plot

5. Create FC Atlas Similarity Matrix

fc_cor_r = lapply(fc_corplot, function(atlas) atlas$sim_mat$r[lower.tri(atlas$sim_mat$r)])
fc_cor_r_cor = lapply(fc_cor_r, function(r1) sapply(fc_cor_r, function(r2) cor.test(r1,r2)))
fc_cor_r_cor_mat = sapply(fc_cor_r_cor, function(atlas) atlas['estimate',])

levelplot(fc_cor_r_cor_mat,col.regions = rev(rainbow(100))[-c(1:20)], 
          region =T, ylab.right = "Pearson correlation", 
          main=list(label=paste('FC Atlas Similarity Matrix')), xlab="",ylab="", 
          scales=list(x=list(at = 1:length(atlases), labels=atlases,rot=90, tck = 0),
                              y=list(at = 1:length(atlases),labels=atlases, tck = 0)))

6. Compute Whole Brain FC-GPS Similarity Matrix Correlations

wbFC_gps_cor = unlist(sapply(fc_cor_r, function(atlas) cor.test(atlas, gps_corplot$r[lower.tri(gps_corplot$r)]))[c('estimate'),])

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls
9: In readChar(file, size, TRUE) : truncating string with embedded nuls

wbFC_gps_cor_df = data_frame(atlas = gsub("\\..*","",names(wbFC_gps_cor)), corr = wbFC_gps_cor)
p = ggplot(data=wbFC_gps_cor_df, aes(x=atlas, y=corr)) +
  geom_bar(stat="identity") + scale_y_continuous(limits=c(-0.5,0.5)) +
  theme_cowplot() + theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) +
  ylab("FC-GPS Similarity Matrices Correlatons") + xlab("")
p

7. Compute Community-Level FC-GPS Similarity Matrix Correlations

There were 13 warnings (use warnings() to see them)

fc_com_gps = function(subj_net_files_1,subj_net_files_2,atlas_now){
  atlas_path = "/Users/hxia/Documents/GitHub/xcpEngine/atlas"
  fc_com_1 = get_fc_com(atlas_now, subj_net = subj_net_files_1)
  fc_com_2 = get_fc_com(atlas_now, subj_net = subj_net_files_2)
  fc_com_corplot = sim_plots(subj_seq = gps_corplot$subj, ids,names(fc_com_1) , fc_com_1, fc_com_2)
  fc_com_cor_r = lapply(fc_com_corplot, function(atlas) atlas$sim_mat$r[lower.tri(atlas$sim_mat$r)])
  fc_com_gps_pvalue = unlist(sapply(fc_com_cor_r, function(atlas) suppressWarnings(cor.test(atlas, gps_corplot$r[lower.tri(gps_corplot$r)], method = "spearman")$p.value)))
  fc_com_gps_cor = unlist(sapply(fc_com_cor_r, function(atlas) suppressWarnings(cor.test(atlas, gps_corplot$r[lower.tri(gps_corplot$r)], method = "spearman")$estimate)))
  fc_com_gps_p_perm = unlist(sapply(fc_com_cor_r, function(atlas) cor_perm(atlas, gps_corplot$r[lower.tri(gps_corplot$r)],1000,  method = "spearman")))
                          
                          
  fc_com_gps_cor_df = data_frame(atlas = gsub("\\..*","",names(fc_com_gps_cor)), corr = fc_com_gps_cor, p_val = fc_com_gps_pvalue, p_perm = fc_com_gps_p_perm )
  fc_com_gps_cor_df$p_val_adj = p.adjust(fc_com_gps_cor_df$p_val,method = "bonferroni")
  fc_com_gps_cor_df$p_perm_adj = p.adjust(fc_com_gps_cor_df$p_perm,method = "bonferroni")
  fc_com_gps_cor_df$sig = fc_com_gps_cor_df$p_val_adj< 0.05
  fc_com_gps_cor_df$sig_perm = fc_com_gps_cor_df$p_perm_adj< 0.05
  q = ggplot(data=fc_com_gps_cor_df, aes(x=reorder(atlas, -corr), y=corr, fill = sig)) +
    geom_bar(stat="identity") + scale_y_continuous(limits=c(-0.5,0.5)) +
    theme_cowplot() + theme(axis.text.x = element_text(angle = 90)) +
    ylab("FC-GPS Similarity Matrices Correlatons") + xlab("") + ggtitle(atlas_now) + theme(plot.title = element_text(hjust = 0.5))
  return(list(fc_corplot = fc_com_corplot, fc_gps_df = fc_com_gps_cor_df, plot = p, com_names = com_names))

fc_schaef400_7_com_gps = fc_com_gps(subj_net_files_1,subj_net_files_2,"schaefer400x7")

[1] "atlas: vis-vis"
[1] 82
[1] "atlas: vis-som"
[1] 82
[1] "atlas: vis-dor"
[1] 82
[1] "atlas: vis-sal"
[1] 82
[1] "atlas: vis-lim"
[1] 82
[1] "atlas: vis-fro"
[1] 82
[1] "atlas: vis-def"
[1] 82
[1] "atlas: som-som"
[1] 82
[1] "atlas: som-dor"
[1] 82
[1] "atlas: som-sal"
[1] 82
[1] "atlas: som-lim"
[1] 82
[1] "atlas: som-fro"
[1] 82
[1] "atlas: som-def"
[1] 82
[1] "atlas: dor-dor"
[1] 82
[1] "atlas: dor-sal"
[1] 82
[1] "atlas: dor-lim"
[1] 82
[1] "atlas: dor-fro"
[1] 82
[1] "atlas: dor-def"
[1] 82
[1] "atlas: sal-sal"
[1] 82
[1] "atlas: sal-lim"
[1] 82
[1] "atlas: sal-fro"
[1] 82
[1] "atlas: sal-def"
[1] 82
[1] "atlas: lim-lim"
[1] 82
[1] "atlas: lim-fro"
[1] 82
[1] "atlas: lim-def"
[1] 82
[1] "atlas: fro-fro"
[1] 82
[1] "atlas: fro-def"
[1] 82
[1] "atlas: def-def"
[1] 82
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000
[1] 1000

fc_schaef200_7_com_gps$plot / fc_schaef400_7_com_gps$plot / fc_power_com_gps$plot

fc_com_gps_out_corplot = function(fc_com_gps_out, title) {
  vector_to_mat(fc_com_gps_out$fc_gps_df$atlas,"half+d")
  com_gps_corr_mat = vector_to_mat(fc_com_gps_out$fc_gps_df$corr,"half+d")
  com_names = fc_com_gps_out$com_names
  rownames(com_gps_corr_mat) = com_names
  colnames(com_gps_corr_mat) = com_names
  com_gps_pval_mat = vector_to_mat(fc_com_gps_out$fc_gps_df$p_val_adj,"half+d")
  col2 <- colorRampPalette(c("#67001F", "#B2182B", 
                             "#FFFFFF", 
                             "#2166AC", "#053061"))
corrplot(com_gps_corr_mat, type="full",method="color", cl.lim = c(-0.2,0.2), tl.col="black", 
           p.mat = com_gps_pval_mat, sig.level = 0.05, insig = "blank", 
           col = rev(col2(200)), title = title)

}

fc_com_gps_out_corplot(fc_power_com_gps, "power264")

Error in dimnames(x) <- dn : 
  length of 'dimnames' [1] not equal to array extent

7. Match FC target to database

8. Permutation Test for Matching

subj_net_perm_1 = get_subj_net_perm(subj_net_files_1)

Error in sample(length(subj_net_files[[atalas]])) : 
  object 'atalas' not found

data_long <- gather(perm_acc, key = "atlas", value = "accuracy", atlases)

Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(atlases)` instead of `atlases` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.

get_com_fc_accuracy = function(atlas_now){
  com_nums = read.table(file.path(atlas_path,atlas_now,paste0(atlas_now,"CommunityAffiliation.1D")))$V1
  com_names = read.table(file.path(atlas_path,atlas_now,paste0(atlas_now,"CommunityNames.txt")))$V1
  com_names = substr(com_names,0,3)
  if (atlas_now == "power264") {com_names[c(1,2)] = c("smH","smM")}
  fc_com_1  = get_fc_com(com_nums, com_names, atlas_now, subj_net_files_1)
  fc_com_2  = get_fc_com(com_nums, com_names, atlas_now, subj_net_files_2)
  #fc_com_1_sig = fc_com_1[names(which(fc_schaef400_7_com_gps$fc_gps_df$sig == T))]
  #fc_com_2_sig = fc_com_2[names(which(fc_schaef400_7_com_gps$fc_gps_df$sig == T))]
  single_multi_acc_com = get_acc_finger(fc_com_1,fc_com_2)
  single_multi_acc_com_df = data.frame(atlas = names(single_multi_acc_com$acc), accuracy = single_multi_acc_com$acc, FC_GPS_sig = fc_schaef400_7_com_gps$fc_gps_df$sig == T)
  p = ggplot(data=single_multi_acc_com_df, aes(x=reorder(atlas, accuracy), y=accuracy, fill = FC_GPS_sig)) +
    geom_bar(stat="identity") +
    theme_cowplot() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
    ylab("FC Fingerprint Accuracy") + xlab("")

return(list(df = single_multi_acc_com_df, plot = p))
}

s400_fc_com_acc$plot / s200_fc_com_acc$plot

atlas_now

[1] "schaefer400x7"

 

---

A work by Cedric Huchuan Xia

hxia@upenn.edu